*** GOF test *******************************************************************

use "$pathD/all_tract.dta", clear

// Generate the temporal fixed effects.
gen month = month(date)
gen dow = dow(date)

// Set the globals.
global weather_covs = "AIR_temp_* AIR_dew_* PRECIP_quant* WIND_speed_*"
global time = "i.month i.dow"
global opt = ", fe cluster(tract_id)"

// Set the panel variable.
destring tract_id, replace
xtset tract_id

// Store the ll estimates.
putexcel set "$pathR/Data/boxcox.xlsx", replace
local j = 1
forvalues i = -1(.1)1 {
	qui: gen indoor_t = (PM_2p5_indoor^`i'-1)/`i'
	qui: xtreg indoor_t PM_2p5_outdoor $weather_covs $time $opt
	qui: drop *_t
	qui: estat ic
	putexcel B`j' = matrix(e(ll))
	local j = `j' + 1
	}
qui: xtreg log_PM_2p5_indoor PM_2p5_outdoor $weather_covs $time $opt
qui: estat ic
putexcel B11 = matrix(e(ll))

// Generate the figure.
import excel "$pathR/Data/boxcox.xlsx", clear
rename B ll_score
gen theta = -1 if _n == 1
replace theta = theta[_n-1] + .1 if _n != 1
twoway (scatter ll_score theta, mcolor(dknavy%100) msize(vsmall) sort) ///
	, ///
	yscale(range(-800000 0) lcolor(gs3)) ///
	ylabel(-800000 "-800K" -600000 "-600K" -400000 "-400K" -200000 "-200K" ///
		0, angle(horizontal) labsize(small) labcolor(gs3) ///
		tlcolor(gs3) glcolor(white)) ///
	ytitle("") ///
	xscale(range(-1 1) lcolor(gs3)) ///
	xlabel(-1(.2)1, nogrid labsize(small)) ///
	xticks(-1(.1)1) ///
	xtitle("") ///
	legend(off) ///
	graphregion(margin(0 0 5 0) fcolor(white) ///
		lcolor(white) ifcolor(white) ilcolor(white)) ///
	text(0 -1 "Log-likelihood score", ///
		color(gs3) size(medsmall) j(left) placement(r)) ///
	text(-1000000 1 "{&theta}", ///
		color(gs3) size(medsmall) j(right) placement(l))
graph export "$pathR/Figures/ll_score.png", replace as(png)

*** Plots **********************************************************************

// Load the data.
use "$pathD/all_tract.dta", clear

// Generate the temporal fixed effects.
gen month = month(date)
gen dow = dow(date)

// Set the globals.
global weather_covs = ///
	"AIR_temp_* AIR_dew_* PRECIP_quant* WIND_speed_*"
global time = "i.month i.dow"
global opt = ", fe cluster(tract_id)"

// Set the panel variable.
destring tract_id, replace
xtset tract_id

// Lin-lin pollution test
preserve
xtreg PM_2p5_indoor PM_2p5_outdoor $weather_covs $time $opt
xtavplot PM_2p5_outdoor, generate(exvar eyvar) nod
xtile pc_exvar = exvar, n(10000)
bysort pc_exvar: egen mean_eyvar = mean(eyvar)
bysort pc_exvar: egen mean_exvar = mean(exvar)
twoway (scatter mean_eyvar mean_exvar, mcolor(dknavy%100) msize(tiny)) ///
	(lfit eyvar exvar, lcolor(cranberry)) ///
	, ///
	yscale(range() lcolor(gs3)) ///
	ylabel(-20(20)80, angle(horizontal) labsize(small) labcolor(gs3) ///
		tlcolor(gs3) glcolor(white)) ///
	ytitle("") ///
	xscale(range() lcolor(gs3)) ///
	xlabel(-100(100)400, labsize(small) labcolor(gs3) ///
		tlcolor(gs3) glcolor(white)) ///
	xtitle("") ///
	legend(off) ///
	graphregion(margin(0 0 5 0) fcolor(white) ///
		lcolor(white) ifcolor(white) ilcolor(white)) ///
	text(80 -100 ///
		"Residual indoor PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(medsmall) j(left) placement(r)) ///
	text(-32 400 ///
		"Residual outdoor PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(medsmall) j(right) placement(l))
graph export "$pathR/Figures/lin_lin_adjusted.png", replace as(png)
restore

// Log-log pollution test
preserve
xtreg log_PM_2p5_indoor log_PM_2p5_outdoor $weather_covs $time $opt
xtavplot log_PM_2p5_outdoor, generate(exvar eyvar) nod
xtile pc_exvar = exvar, n(10000)
bysort pc_exvar: egen mean_eyvar = mean(eyvar)
bysort pc_exvar: egen mean_exvar = mean(exvar)
twoway (scatter mean_eyvar mean_exvar, mcolor(dknavy%100) msize(tiny)) ///
	(lfit eyvar exvar, lcolor(cranberry)) ///
	, ///
	yscale(range() lcolor(gs3)) ///
	ylabel(-6(2)4, angle(horizontal) labsize(small) labcolor(gs3) ///
		tlcolor(gs3) glcolor(white)) ///
	ytitle("") ///
	xscale(range() lcolor(gs3)) ///
	xlabel(-8(2)6, labsize(small) labcolor(gs3) ///
		tlcolor(gs3) glcolor(white)) ///
	xtitle("") ///
	legend(off) ///
	graphregion(margin(0 0 5 0) fcolor(white) ///
		lcolor(white) ifcolor(white) ilcolor(white)) ///
	text(4 -8 ///
		"Residual log indoor PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(medsmall) j(left) placement(r)) ///
	text(-7.2 6 ///
		"Residual log outdoor PM{subscript:2.5} ({&mu}g/m{superscript:3})", ///
		color(gs3) size(medsmall) j(right) placement(l))
graph export "$pathR/Figures/log_log_adjusted.png", replace as(png)
restore

clear
